#ifndef __CURLTP13D3D__
#define __CURLTP13D3D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// Hcurl P1 element, CURL CONFORMING EDGE ELEMENTS OF NÉDÉLEC, 3D
// ***********************************************************************
static INT Curl_T_P1_3D_3D_dof[4] = {0, 1, 0, 0};
static INT Curl_T_P1_3D_3D_Num_Bas = 6;
static INT Curl_T_P1_3D_3D_Value_Dim = 3;
static INT Curl_T_P1_3D_3D_Inter_Dim = 1;
static INT Curl_T_P1_3D_3D_Polydeg = 1;
static bool Curl_T_P1_3D_3D_IsOnlyDependOnRefCoord = 0;
static INT Curl_T_P1_3D_3D_Accuracy = 1;
static MAPTYPE Curl_T_P1_3D_3D_Maptype = Piola;

static void Curl_T_P1_3D_3D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 1 - z - y;
    values[1] = x;
    values[2] = x;
    values[3] = y;
    values[4] = 1 - z - x;
    values[5] = y;
    values[6] = z;
    values[7] = z;
    values[8] = 1 - y - x;
    values[9] = -y;
    values[10] = x;
    values[11] = 0;
    values[12] = -z;
    values[13] = 0;
    values[14] = x;
    values[15] = 0;
    values[16] = -z;
    values[17] = y;
}
static void Curl_T_P1_3D_3D_D000(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 1 - z - y;
    values[1] = x;
    values[2] = x;
    values[3] = y;
    values[4] = 1 - z - x;
    values[5] = y;
    values[6] = z;
    values[7] = z;
    values[8] = 1 - y - x;
    values[9] = -y;
    values[10] = x;
    values[11] = 0;
    values[12] = -z;
    values[13] = 0;
    values[14] = x;
    values[15] = 0;
    values[16] = -z;
    values[17] = y;
}
static void Curl_T_P1_3D_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 1;
    values[2] = 1;
    values[3] = 0;
    values[4] = -1;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = -1;
    values[9] = 0;
    values[10] = 1;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 1;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = -1;
    values[1] = 0;
    values[2] = 0;
    values[3] = 1;
    values[4] = 0;
    values[5] = 1;
    values[6] = 0;
    values[7] = 0;
    values[8] = -1;
    values[9] = -1;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 1;
}
static void Curl_T_P1_3D_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = -1;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = -1;
    values[5] = 0;
    values[6] = 1;
    values[7] = 1;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = -1;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = -1;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x = RefCoord[0];
    y = RefCoord[1];
    z = RefCoord[2];
    values[0] = 0;
    values[1] = 0;
    values[2] = 0;
    values[3] = 0;
    values[4] = 0;
    values[5] = 0;
    values[6] = 0;
    values[7] = 0;
    values[8] = 0;
    values[9] = 0;
    values[10] = 0;
    values[11] = 0;
    values[12] = 0;
    values[13] = 0;
    values[14] = 0;
    values[15] = 0;
    values[16] = 0;
    values[17] = 0;
}
static void Curl_T_P1_3D_3D_Nodal(ELEMENT *elem, FUNCTIONVEC *fun, INT dim, DOUBLE *values)
{
    // 直接给出积分格式
    INT NumPoints = 4;
    DOUBLE QuadX14[4] = {0.069431844202973713731097404888715,
                         0.33000947820757187134432797392947,
                         0.66999052179242812865567202607053,
                         0.93056815579702628626890259511129};
    DOUBLE QuadW14[4] = {0.17392742256872692485636378,
                         0.3260725774312730473880606,
                         0.3260725774312730473880606,
                         0.17392742256872692485636378};
    INT vert0[6] = {0, 0, 0, 1, 1, 2};
    INT vert1[6] = {1, 2, 3, 2, 3, 3};
    INT i, j, k;
    INT num_fun = dim / 3;
    DOUBLE coord[3];
    DOUBLE vec_line[3]; // 带方向的线(按照单元内部的方向)
    DOUBLE *values_temp = malloc(dim * sizeof(DOUBLE));
#if NEDELECSCALE == 0
    for (i = 0; i < 6; i++)
    {
        for (k = 0; k < num_fun; k++)
            values[k + i * num_fun] = 0.0;
        vec_line[0] = elem->Vert_X[vert1[i]] - elem->Vert_X[vert0[i]];
        vec_line[1] = elem->Vert_Y[vert1[i]] - elem->Vert_Y[vert0[i]];
        vec_line[2] = elem->Vert_Z[vert1[i]] - elem->Vert_Z[vert0[i]];
        for (j = 0; j < NumPoints; j++)
        {
            // coord
            coord[0] = QuadX14[j] * elem->Vert_X[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_X[vert1[i]];
            coord[1] = QuadX14[j] * elem->Vert_Y[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_Y[vert1[i]];
            coord[2] = QuadX14[j] * elem->Vert_Z[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_Z[vert1[i]];
            // u(t_j)
            fun(coord, dim, values_temp);
            // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau}\omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j
            for (k = 0; k < num_fun; k++)
            {
                values[k + i * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j];
            }
        }
    }
#else
    DOUBLE length;
    for (i = 0; i < 6; i++)
    {
        for (k = 0; k < num_fun; k++)
            values[k + i * num_fun] = 0.0;
        vec_line[0] = elem->Vert_X[vert1[i]] - elem->Vert_X[vert0[i]];
        vec_line[1] = elem->Vert_Y[vert1[i]] - elem->Vert_Y[vert0[i]];
        vec_line[2] = elem->Vert_Z[vert1[i]] - elem->Vert_Z[vert0[i]];
        length = sqrt(vec_line[0]*vec_line[0]+vec_line[1]*vec_line[1]+vec_line[2]*vec_line[2]);
        for (j = 0; j < NumPoints; j++)
        {
            // coord
            coord[0] = QuadX14[j] * elem->Vert_X[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_X[vert1[i]];
            coord[1] = QuadX14[j] * elem->Vert_Y[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_Y[vert1[i]];
            coord[2] = QuadX14[j] * elem->Vert_Z[vert0[i]] + (1 - QuadX14[j]) * elem->Vert_Z[vert1[i]];
            // u(t_j)
            fun(coord, dim, values_temp);
            // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau}\omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j
            for (k = 0; k < num_fun; k++)
            {
                values[k + i * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j] / length;
            }
        }
    }
#endif
    free(values_temp);
}
#endif
